library(prism)
library(raster)
library(lubridate)
library(geosphere)
library(RColorBrewer)
library(png)
library(grid)
library(rasterVis)
library(rgdal)
library(scales)
library(viridis)
library(ggthemes)
library(tidyverse)
library(lubridate)
library(DomoR)
library(choroplethr)
#get_prism_dailys(type="tmax", minDate = "2017-01-01", maxDate = "2018-11-29", keepZip=F)
#get_prism_dailys(type="tmin", minDate = "2017-01-01", maxDate = "2018-11-29", keepZip=F)
#get_prism_dailys(type="tmean", minDate = "2017-01-01", maxDate = "2018-11-29", keepZip=F)
#get_prism_dailys(type="ppt", minDate = "2017-01-01", maxDate = "2018-11-29", keepZip=F)
retailCal <- read.csv("C:/Users/sewing/OneDrive - Do My Own/NOAA/4-4-5-calendar.csv",
stringsAsFactors = FALSE) %>%
mutate(Actual.Date = date(mdy(Actual.Date)),
X2018.equiv = date(mdy(X2018.equiv))) %>%
filter(Retail.Year > 2016,
Retail.Year < 2019) %>%
select(Actual.Date, Retail.Week, Retail.Year, Retail.DayofYear) %>%
spread(Retail.Year, Actual.Date) %>%
mutate(meanFile2017 = character(1),
meanFile2018 = character(1),
meanPath2017 = character(1),
meanPath2018 = character(1)) %>%
filter(!is.na(`2017`))
RetailWeeks2018 <- read.csv("C:/Users/sewing/OneDrive - Do My Own/NOAA/4-4-5-calendar.csv",
stringsAsFactors = FALSE) %>%
filter(Retail.Year == 2018) %>%
select(Retail.Week, Actual.Date) %>%
mutate(Actual.Date = mdy(Actual.Date)) %>%
group_by(Retail.Week) %>%
summarise(max_date = max(Actual.Date),
min_date = min(Actual.Date)) %>%
mutate(min_month = as.character(month(min_date, label = TRUE)),
max_month = as.character(month(max_date, label = TRUE)),
months = ifelse(min_month == max_month, min_month, paste0(min_month, "/", max_month))) %>%
select(Retail.Week, months)
options(prism.path = "C:/Users/sewing/OneDrive - Do My Own/NOAA/PRISM_Files")
prismFileList <- ls_prism_data(absPath = TRUE)
stableMeans <- prismFileList[which(grepl("_tmean_stable_4kmD1_",
prismFileList$files)),] %>%
mutate(Date = ymd(str_sub(files, start = 26, end = -5)))
lastStableDate <- max(stableMeans$Date)
get_prism_dailys(type="tmean",
minDate = lastStableDate,
maxDate = Sys.Date(),
keepZip=F)
prismFileList <- ls_prism_data(absPath = TRUE)
stableMeans <- prismFileList[which(grepl("_tmean_stable_4kmD1_",
prismFileList$files)),] %>%
mutate(Date = ymd(str_sub(files, start = 26, end = -5)))
newLastStableDate <- max(stableMeans$Date)
provisionalMeans <- prismFileList[which(grepl("_tmean_provisional_4kmD1_",
prismFileList$files)),] %>%
mutate(Date = ymd(str_sub(files, start = -12, end = -5))) %>%
filter(Date > ymd(newLastStableDate))
lastProvisionalDate <- max(provisionalMeans$Date)
earlyMeans <- prismFileList[which(grepl("_tmean_early_4kmD1_",
prismFileList$files)),] %>%
mutate(Date = ymd(str_sub(files, start = -12, end = -5))) %>%
filter(Date > ymd(lastProvisionalDate))
lastEarlyDate <- max(earlyMeans$Date)
meansFiles <- stableMeans %>%
rbind(provisionalMeans) %>%
rbind(earlyMeans)
prismComps <- retailCal %>%
mutate(meanFile2017 = meansFiles$files[which(meansFiles$Date == retailCal$`2017`)],
meanFile2018 = ifelse(`2018` > lastEarlyDate,
"Future Date",
meansFiles$files[which(meansFiles$Date == retailCal$`2018`)]),
meanPath2017 = meansFiles$abs_path[which(meansFiles$Date == retailCal$`2017`)],
meanPath2018 = ifelse(`2018` > lastEarlyDate,
"Future Date",
meansFiles$abs_path[which(meansFiles$Date == retailCal$`2018`)]))
dailyRasters <- prismComps %>%
mutate(`2017 Rasters` = list(length = 1),
`2018 Rasters` = list(length = 1))
i <- 1
for(i in 1:dim(dailyRasters)[1]){
dailyRasters$`2017 Rasters`[i] <- raster(dailyRasters$meanPath2017[i])
i = i + 1
}
i <- 1
for(i in 1:table(dailyRasters$`2018`<lastEarlyDate)[2]){
dailyRasters$`2018 Rasters`[i] <- raster(dailyRasters$meanPath2018[i])
i = i + 1
}
anomCalc <- function(x, y) {
return(x - y)
}
avgWeek <- function(a, b, c, d, e, f, g) {
return((((a+b+c+d+e+f+g)/7)*(9/5))+32)
}
i <- 0
while(i < week(lastEarlyDate)-1){
i = i+1
weekOverlay17 <- overlay(dailyRasters$`2017 Rasters`[[(7*(i-1))+1]],
dailyRasters$`2017 Rasters`[[(7*(i-1))+2]],
dailyRasters$`2017 Rasters`[[(7*(i-1))+3]],
dailyRasters$`2017 Rasters`[[(7*(i-1))+4]],
dailyRasters$`2017 Rasters`[[(7*(i-1))+5]],
dailyRasters$`2017 Rasters`[[(7*(i-1))+6]],
dailyRasters$`2017 Rasters`[[(7*(i-1))+7]],
fun = avgWeek)
writeRaster(weekOverlay17,
filename = paste0("C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/Week_", i, "_2017_avg_temp"),
overwrite=TRUE)
weekOverlay18 <- overlay(dailyRasters$`2018 Rasters`[[(7*(i-1))+1]],
dailyRasters$`2018 Rasters`[[(7*(i-1))+2]],
dailyRasters$`2018 Rasters`[[(7*(i-1))+3]],
dailyRasters$`2018 Rasters`[[(7*(i-1))+4]],
dailyRasters$`2018 Rasters`[[(7*(i-1))+5]],
dailyRasters$`2018 Rasters`[[(7*(i-1))+6]],
dailyRasters$`2018 Rasters`[[(7*(i-1))+7]],
fun = avgWeek)
writeRaster(weekOverlay18,
filename = paste0("C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/Week_", i, "_2018_avg_temp"),
overwrite=TRUE)
weekYoYDiff <- overlay(weekOverlay18, weekOverlay17, fun = anomCalc)
writeRaster(weekYoYDiff,
filename = paste0("C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/Week_", i, "_YoY_diff_avg_temp"),
overwrite=TRUE)
}
YoYFiles <- dir(path = "C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/",
pattern = "YoY_diff_avg_temp.grd",
include.dirs = TRUE,
full.names = TRUE)
us_map <- map_data("state")
i <- 1
while(i < week(Sys.Date())){
r <- raster(paste0("C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Rasters/Week_",
i,
"_YoY_diff_avg_temp.grd"))
r_spdf <- as(r, "SpatialPixelsDataFrame")
r_df <- as.data.frame(r_spdf)
colnames(r_df) <- c("value", "x", "y")
rp <- ggplot() +
geom_path(
data = us_map,
aes(x = long,
y = lat,
group = group),
color = "white",
size = 0.1) +
coord_equal() +
theme_map() +
theme(legend.position = "bottom") +
theme(legend.key.width = unit(2, "cm")) +
labs(
x = NULL,
y = NULL,
title = paste0(RetailWeeks2018[i, 2], " - Week ", i, " YoY Diff in Temperature")
) +
geom_raster(data = r_df, aes(x = x, y = y, fill = value)) +
scale_fill_distiller(palette = "RdYlBu", limits = c(-30, 30))
fname <- paste0('C:/Users/sewing/OneDrive - Do My Own/NOAA/Weekly_Plots/Week_',
i,
'_2018_less_2017_Plot.png')
ggsave(rp, file=fname)
i = i+1
}